rm(list = ls()) # Borramos variables de environment
# descomentar si es la primera vez (y requieren instalación)
# install.packages("tidyverse")
# install.packages("performance")
# install.packages("olsrr")
# install.packages("corrr")
# install.packages("corrplot")
# install.packages("skimr")
library(tidyverse)
library(rsample)
library(performance)
library(olsrr)
library(corrr)
library(corrplot)
library(skimr)Entrega II
Instrucciones (leer antes de empezar)
Modifica de la cabecera del documento
.qmdtus datos personales (nombre y DNI). IMPORTANTE: no modifiques nada más de la cabecera.Asegúrate, ANTES de seguir editando el documento, que el archivo
.qmdse renderiza correctamente y se genera el.htmlcorrespondiente en tu carpeta local de tu ordenador (pulsa el botónRendero guarda el documento conRender on Saveactivado).Los chunks (cajas de código) creados están o vacíos o incompletos, de ahí que la mayoría tengan la opción
#| eval: false. Una vez que edites lo que consideres, debes ir cambiando cada chunck a#| eval: true(o quitarlo directamente) para que se ejecuten.Recuerda que puedes ejecutar chunk a chunk con el botón play o ejecutar todos los chunk hasta uno dado (con el botón a la izquierda del anterior).
IMPORTANTE: solo se podrá subir un archivo
.htmlal campus, no se evaluarán entregas sin el.htmlcorrectamente generado. Recuerda: el código es una herramienta, no el fin en esta asignatura. Se evaluará especialmente las interpretaciones y conclusiones que detalles. Un ejercicio con el código perfecto pero sin ningún tipo de razonamiento, interpretación o conclusión no superará el 4 sobre 10 de nota.
Paquetes necesarios
Necesitaremos los siguientes paquetes (haz play en el chunk para que se carguen):
Caso práctico: análisis de datos de cáncer
El archivo de datos a usar lo cargaremos desde el csv cancer.csv.
# no cambies código
datos <- read_csv(file = "./cancer.csv")
datos# A tibble: 3,047 × 13
deathRate medIncome popEst2015 povertyPercent studyPerCap MedianAge
<dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 165. 61898 260131 11.2 500. 39.3
2 161. 48127 43269 18.6 23.1 33
3 175. 49348 21026 14.6 47.6 45
4 195. 44243 75882 17.1 343. 42.8
5 144. 49955 10321 12.5 0 48.3
6 176 52313 61023 15.6 180. 45.4
7 176. 37782 41516 23.2 0 42.6
8 184. 40189 20848 17.8 0 51.7
9 190. 42579 13088 22.3 0 49.3
10 178. 60397 843954 13.1 428. 35.8
# ℹ 3,037 more rows
# ℹ 7 more variables: AvgHouseholdSize <dbl>, PercentMarried <dbl>,
# PctHS25_Over <dbl>, PctUnemployed16_Over <dbl>, PctPrivateCoverage <dbl>,
# BirthRate <dbl>, region <chr>Los datos representan diferentes características socioeconómicas de distintas regiones, extraídas de the American Community Survey (census.gov), clinicaltrials.gov, y cancer.gov.
¿Nuestro objetivo? Ser capaces de predecir de manera lineal y MULTIVARIANTE la variable deathRate, nuestra variable objetivo, que representa la mortalidad media de cancer, por cada 100 000 habitantes. El resto de variables son:
medianIncome: mediana de los ingresos de la región.popEst2015: población de la regiónpovertyPercent: porcentaje de población en situación de pobreza.studyPerCap: ensayos clínicos relacionados por el cáncer realizados por cada 100 000 habitantes.MedianAge: edad mediana.region: nombre de la región.AvgHouseholdSize: tamaño medio de los hogares.PercentMarried: porcentaje de habitantes casados.PctHS25_Over: porcentaje de residentes por encima de los 25 años con (máximo) título de bachillerato.PctUnemployed16_Over: porcentaje de residentes de más de 16 años en paro.PctPrivateCoverage: porcentaje de residentes con seguro de salud privado.BirthRate: tasa de natalidad.
IMPORTANTE: siempre que puedas, relaciona los valores numéricos de la salida de R con su fórmula.
Ejercicio 1 (1.25 puntos)
Chequea que no hay problemas de rango/codificación: en caso de que alguna variable tenga dichos problemas, debes sustituir los valores fueran del rango razonable por su media o mediana (sin esos datos fuera del rango), justificando por qué usas la media o por qué usas la mediana, ejecutando el código que consideres (si te atascas, sigue haciendo con el dataset original)
Como se vio en la entrega anterior, por los rangos aportados, el único problema de rango/codificación lo tenemos en MedianAge, con valores por encima de un rango razonable.
datos |> skim()| Name | datos |
| Number of rows | 3047 |
| Number of columns | 13 |
| _______________________ | |
| Column type frequency: | |
| character | 1 |
| numeric | 12 |
| ________________________ | |
| Group variables | None |
Variable type: character
| skim_variable | n_missing | complete_rate | min | max | empty | n_unique | whitespace |
|---|---|---|---|---|---|---|---|
| region | 0 | 1 | 4 | 9 | 0 | 4 | 0 |
Variable type: numeric
| skim_variable | n_missing | complete_rate | mean | sd | p0 | p25 | p50 | p75 | p100 | hist |
|---|---|---|---|---|---|---|---|---|---|---|
| deathRate | 0 | 1 | 178.66 | 27.75 | 59.70 | 161.20 | 178.10 | 195.20 | 362.80 | ▁▇▆▁▁ |
| medIncome | 0 | 1 | 47063.28 | 12040.09 | 22640.00 | 38882.50 | 45207.00 | 52492.00 | 125635.00 | ▇▇▁▁▁ |
| popEst2015 | 0 | 1 | 102637.37 | 329059.22 | 827.00 | 11684.00 | 26643.00 | 68671.00 | 10170292.00 | ▇▁▁▁▁ |
| povertyPercent | 0 | 1 | 16.88 | 6.41 | 3.20 | 12.15 | 15.90 | 20.40 | 47.40 | ▃▇▃▁▁ |
| studyPerCap | 0 | 1 | 155.40 | 529.63 | 0.00 | 0.00 | 0.00 | 83.65 | 9762.31 | ▇▁▁▁▁ |
| MedianAge | 0 | 1 | 45.27 | 45.30 | 22.30 | 37.70 | 41.00 | 44.00 | 624.00 | ▇▁▁▁▁ |
| AvgHouseholdSize | 0 | 1 | 2.48 | 0.43 | 0.02 | 2.37 | 2.50 | 2.63 | 3.97 | ▁▁▃▇▁ |
| PercentMarried | 0 | 1 | 51.77 | 6.90 | 23.10 | 47.75 | 52.40 | 56.40 | 72.50 | ▁▂▇▇▁ |
| PctHS25_Over | 0 | 1 | 34.80 | 7.03 | 7.50 | 30.40 | 35.30 | 39.65 | 54.80 | ▁▂▇▇▁ |
| PctUnemployed16_Over | 0 | 1 | 7.85 | 3.45 | 0.40 | 5.50 | 7.60 | 9.70 | 29.40 | ▅▇▁▁▁ |
| PctPrivateCoverage | 0 | 1 | 64.35 | 10.65 | 22.30 | 57.20 | 65.10 | 72.10 | 92.30 | ▁▂▇▇▂ |
| BirthRate | 0 | 1 | 5.64 | 1.99 | 0.00 | 4.52 | 5.38 | 6.49 | 21.33 | ▂▇▁▁▁ |
ggplot(datos) +
geom_histogram(aes(x = MedianAge)) +
theme_minimal()`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
En el histograma vemos que hay un corte muy claro, por lo que estableceremos (por ejemplo) un máximo de 120 años: por encima de dicha cantidad, se considerará una observación mal codificada. En lugar de eliminarlas lo que vamos a hacer es sustituirla por una medida de centralidad (imputación), que podemos hacer en base a dos argumentos:
- Si usamos todos los datos, el histograma anterior, tenemos una distribución muy asimétrica, por lo que deberíamos usar la mediana (que se ve menos afectada por esos valores mal codificados, ya que es más robusta)
ggplot(datos) +
geom_density(aes(x = MedianAge)) +
theme_minimal()
ggplot(datos) +
geom_boxplot(aes(x = MedianAge)) +
theme_minimal()
- En esta solución usaremos la media PERO sin esos datos (ya que en realidad…deberían ser
NAya que están mal recopilados).
ggplot(datos |> filter(MedianAge < 120)) +
geom_histogram(aes(x = MedianAge)) +
theme_minimal()`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
ggplot(datos |> filter(MedianAge < 120)) +
geom_density(aes(x = MedianAge)) +
theme_minimal()
ggplot(datos |> filter(MedianAge < 120)) +
geom_boxplot(aes(x = MedianAge)) +
theme_minimal()
Vemos en los gráficos que es bastante simétrica y, de hecho, su distribución se asemeja mucho a una normal, por lo que la media (sin esos datos) si es representativa de la población, y lo que haremos será sustituir los valores por encima del límite permitido por dicha media
media <-
datos |>
filter(MedianAge < 120) |>
summarise(mean = mean(MedianAge)) |>
pull(mean)
# Escribe el código que consideres
datos_preproc <-
datos |>
mutate(MedianAge = if_else(MedianAge < 120, MedianAge,
media))Ejercicio 2 (0.75 puntos)
En caso de encontrar alguna variable cualitativa, procesa los datos como consideras para incluir su información de manera que pueda ser usada en la futura regresión. Recuerda: si tenemos una variable cualitativa que toma k valores, debemos crear k-1 nuevas variables numéricas (de un tipo muy concreto) en su lugar.
Sí tenemos una variable cualitativa que toma 4 categorías (cualitativa nominal)
datos_preproc |> count(region)# A tibble: 4 × 2
region n
<chr> <int>
1 Midwest 1021
2 Northeast 210
3 South 1368
4 West 448Para poder incluirla en la regresión la recodificaremos con la técnica conocida como one-hot encoding: creamos nuevas variables binarias (si tenemos \(k=4\) modalidades, crearemos \(k-1 = 3\) variables) que nos diga si el registro toma dicho valor o no (hacemos una menos para evitar colinealidad)
# Escribe el código que consideres
datos_preproc2 <-
datos_preproc |>
mutate(region_Northeast = region == "Northeast",
region_south = region == "South",
region_West = region == "West") |>
select(-region)Esto lo podemos hacer de manera automática con dummy_cols del paquete {fastDummies}
# Escribe el código que consideres
datos_preproc2 <-
datos_preproc |>
fastDummies::dummy_cols(select_columns = "region",
remove_first_dummy = TRUE) |>
select(-region)Ejercicio 3 (1.5 puntos)
Ejecuta el código que consideres para realizar una selección inicial de variables de manera que se mitiguen los efectos adversos de la colinealidad. Comenta todo lo que consideres sobre. Tras ello separa las observaciones en dos datasets distintos: un dataset train y un segundo dataset test. Para ello realiza un muestreo aleatorio de manera que train contenga el 80% de los datos y test el 20% restante. IMPORTANTE: no cambies la semilla.
Vamos a realizar un análisis de correlaciones previo
datos_preproc2 |>
corrr::correlate()Correlation computed with
• Method: 'pearson'
• Missing treated using: 'pairwise.complete.obs'# A tibble: 15 × 16
term deathRate medIncome popEst2015 povertyPercent studyPerCap MedianAge
<chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 deathRate NA -0.429 -0.120 0.429 -0.0223 -0.00426
2 medIncome -0.429 NA 0.236 -0.789 0.0440 -0.117
3 popEst20… -0.120 0.236 NA -0.0653 0.0557 -0.176
4 povertyP… 0.429 -0.789 -0.0653 NA -0.0557 -0.193
5 studyPer… -0.0223 0.0440 0.0557 -0.0557 NA -0.0317
6 MedianAge -0.00426 -0.117 -0.176 -0.193 -0.0317 NA
7 AvgHouse… -0.0369 0.112 0.110 0.0743 -0.00407 -0.356
8 PercentM… -0.267 0.355 -0.160 -0.643 -0.0381 0.429
9 PctHS25_… 0.405 -0.471 -0.312 0.194 -0.0851 0.330
10 PctUnemp… 0.378 -0.453 0.0508 0.655 -0.0320 -0.127
11 PctPriva… -0.386 0.724 0.0527 -0.823 0.0925 0.0687
12 BirthRate -0.0874 -0.0102 -0.0577 -0.0123 0.0107 -0.0999
13 region_N… -0.0646 0.195 0.135 -0.152 0.0143 0.0621
14 region_S… 0.331 -0.300 -0.0548 0.421 -0.108 -0.118
15 region_W… -0.279 0.114 0.0897 -0.0530 -0.00804 -0.0405
# ℹ 9 more variables: AvgHouseholdSize <dbl>, PercentMarried <dbl>,
# PctHS25_Over <dbl>, PctUnemployed16_Over <dbl>, PctPrivateCoverage <dbl>,
# BirthRate <dbl>, region_Northeast <dbl>, region_South <dbl>,
# region_West <dbl>datos_preproc2 |>
cor() |>
corrplot::corrplot(method = "color")
Concretamente vamos a fijarnos en las correlaciones vs la variable objetivo
datos_preproc2 |>
corrr::correlate() |>
corrr::focus("deathRate") |>
filter(abs(deathRate) < 0.05)Correlation computed with
• Method: 'pearson'
• Missing treated using: 'pairwise.complete.obs'# A tibble: 3 × 2
term deathRate
<chr> <dbl>
1 studyPerCap -0.0223
2 MedianAge -0.00426
3 AvgHouseholdSize -0.0369 Como vemos tenemos 3 variables con una correlación con la variable objetivo por debajo de 0.05 en valor absoluto (prácticamente incorreladas): studyPerCap, MedianAge y AvgHouseholdSize. Podríamos no haberlo hecho (ya que no se pedía el enunciado pero si lo has hecho, es lo correcto) o incluso ser más exigente y subir el umbral de correlación. Vamos a quitar de momento solo esas 3 que presentan una aparente incorrelación lineal (!= independencia) con la objetivo
datos_cor <-
datos_preproc2 |>
select(-studyPerCap, -MedianAge, -AvgHouseholdSize)
datos_cor |>
cor() |>
corrplot::corrplot(method = "color")
Tras ello vemos que hay predictoras altamente correlacionadas entre sí:
medIncomeconpovertyPercentcon $-0.789 $ (lógico: más ingresos, menos pobre)PctPrivateCoverageymedIncomecon \(0.724\)PctPrivateCoverageypovertyPercentcon \(-0.822\)
Para decidir cuál quitamos vamos a basarnos en dos criterios: aquella cuya relación lineal con el conjunto de las demás sea más alta y que tenga menos dependencia con la variable objetivo (tienen más o menos similar entre \(0.43\) y \(0.38\) en valor absoluto)
Este análisis de correlaciones solo nos permite ver relaciones 2 a 2 entre las variables, pero no dependencias lineales más complejas, y para eso vamos a realizar un análisis de multicolinealidad con check_collinearity()
ajuste_saturado <- lm(data = datos_cor, formula = deathRate ~ .)
performance::check_collinearity(ajuste_saturado)# Check for Multicollinearity
Low Correlation
Term VIF VIF 95% CI Increased SE Tolerance Tolerance 95% CI
medIncome 4.22 [3.97, 4.50] 2.06 0.24 [0.22, 0.25]
popEst2015 1.20 [1.15, 1.25] 1.09 0.84 [0.80, 0.87]
PercentMarried 2.30 [2.18, 2.44] 1.52 0.43 [0.41, 0.46]
PctHS25_Over 1.72 [1.64, 1.81] 1.31 0.58 [0.55, 0.61]
PctUnemployed16_Over 2.06 [1.95, 2.17] 1.43 0.49 [0.46, 0.51]
PctPrivateCoverage 4.20 [3.94, 4.47] 2.05 0.24 [0.22, 0.25]
BirthRate 1.07 [1.04, 1.13] 1.04 0.93 [0.89, 0.96]
region_Northeast 1.27 [1.22, 1.33] 1.12 0.79 [0.75, 0.82]
region_South 1.83 [1.74, 1.93] 1.35 0.55 [0.52, 0.58]
region_West 1.60 [1.53, 1.68] 1.27 0.62 [0.59, 0.65]
Moderate Correlation
Term VIF VIF 95% CI Increased SE Tolerance Tolerance 95% CI
povertyPercent 6.83 [6.40, 7.29] 2.61 0.15 [0.14, 0.16]Parece que de lsa 3 la que “colisiona” más con el resto de predictoras es povertyPercent por lo que procederemos retirarla
datos_cor <-
datos_cor |>
select(-povertyPercent)
ajuste_saturado <- lm(data = datos_cor, formula = deathRate ~ .)
performance::check_collinearity(ajuste_saturado)# Check for Multicollinearity
Low Correlation
Term VIF VIF 95% CI Increased SE Tolerance Tolerance 95% CI
medIncome 3.05 [2.88, 3.24] 1.75 0.33 [0.31, 0.35]
popEst2015 1.20 [1.15, 1.25] 1.09 0.84 [0.80, 0.87]
PercentMarried 1.76 [1.67, 1.85] 1.33 0.57 [0.54, 0.60]
PctHS25_Over 1.69 [1.61, 1.78] 1.30 0.59 [0.56, 0.62]
PctUnemployed16_Over 1.98 [1.89, 2.10] 1.41 0.50 [0.48, 0.53]
PctPrivateCoverage 3.48 [3.28, 3.70] 1.87 0.29 [0.27, 0.31]
BirthRate 1.07 [1.04, 1.13] 1.04 0.93 [0.89, 0.96]
region_Northeast 1.26 [1.21, 1.32] 1.12 0.79 [0.76, 0.83]
region_South 1.82 [1.73, 1.91] 1.35 0.55 [0.52, 0.58]
region_West 1.59 [1.52, 1.67] 1.26 0.63 [0.60, 0.66]Por último vamos a realizar la partición pedida
# Escribe el código que consideres
set.seed(12345)
datos_split <- initial_split(datos_cor, prop = 0.8)
train <- training(datos_split)
test <- testing(datos_split)Ejercicio 4 (1 punto)
No usaremos el dataset test hasta el final para evaluar las predicciones (con un dataset que el modelo no conoce). Con el dataset train ejecuta el ajuste de regresión lineal multivariante saturado (habiendo hecho lo pedido en los ejercicios anteriores) y comenta de manera MUY RESUMIDA la salida (SOLO lo que puedas interpretar hasta ahora, y al grano, que se te acaba el tiempo)
# Escribe el código que consideres
ajuste_saturado <- lm(data = train, formula = deathRate ~ .)
ajuste_saturado |> summary()
Call:
lm(formula = deathRate ~ ., data = train)
Residuals:
Min 1Q Median 3Q Max
-103.288 -11.893 0.604 11.863 158.564
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.605e+02 7.315e+00 21.941 < 2e-16 ***
medIncome -1.690e-04 6.480e-05 -2.608 0.00917 **
popEst2015 -1.060e-06 1.352e-06 -0.784 0.43314
PercentMarried -4.777e-01 8.531e-02 -5.599 2.40e-08 ***
PctHS25_Over 1.228e+00 8.342e-02 14.721 < 2e-16 ***
PctUnemployed16_Over 1.669e+00 1.834e-01 9.096 < 2e-16 ***
PctPrivateCoverage -5.346e-02 7.747e-02 -0.690 0.49016
BirthRate -5.073e-01 2.342e-01 -2.167 0.03036 *
region_Northeast -3.588e+00 1.956e+00 -1.835 0.06667 .
region_South 7.577e+00 1.216e+00 6.229 5.51e-10 ***
region_West -1.168e+01 1.602e+00 -7.290 4.18e-13 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 21.99 on 2426 degrees of freedom
Multiple R-squared: 0.3736, Adjusted R-squared: 0.371
F-statistic: 144.7 on 10 and 2426 DF, p-value: < 2.2e-16Algunos comentarios que se podrían hacer:
El \(R^2\) es 0.3736 y el ajustado \(0.371\), lo que nos da un indicio de que tampoco tenemos un sobreajuste muy grande (tras quitar lo que hemos quitado ya, claro)
Ese \(R^2\) prácticamente duplica el univariante de la primera entrega
A mayores ingresos, disminuye la tasa de mortalidad: por cada 10000$ que suba los ingresos medianos, la tasa de mortalidad baja 1.69 puntos.
La educación y empleo importa tanto que por cada % que se incrementa en
PctHS25_Over(es decir, más gente cuyo máximo título es solo bachillerato) o por cada % que se incrementaPctUnemployed16_Over, la tasa de mortalidad sube 1.2 y 1.66 puntos respectivamente (menos formación o menos empleo –> peores condiciones de vida –> más mortalidad)El código postal importa: simplemente por el hecho de estar en el Sur, tu tasa de mortalidad sube ¡7 puntos!
NO SE PUEDE COMENTAR LOS P-VALORES YA QUE NO SABEMOS SI LA DIAGNOSIS ES CORRECTA
Ejercicio 5 (1.75 puntos)
Realiza una selección de modelos en base a los criterios BIC y AIC. Recuerda que para evitar incompatibilidad con otros paquetes, no debes hAcer library(MASS) sino MASS::… Comenta e interpreta todo lo que sepas de las salidas generadas. ¿Cuál penaliza más? ¿Por qué? ¿Qué ventajas tiene la selección de BIC? Interpreta los coeficientes de ambos modelos y comenta las diferencias (en caso de que existiesen)
# Escribe el código que consideres
ajuste_AIC <- MASS::stepAIC(ajuste_saturado, direction = "both",
k = 2)Start: AIC=15073.65
deathRate ~ medIncome + popEst2015 + PercentMarried + PctHS25_Over +
PctUnemployed16_Over + PctPrivateCoverage + BirthRate + region_Northeast +
region_South + region_West
Df Sum of Sq RSS AIC
- PctPrivateCoverage 1 230 1172938 15072
- popEst2015 1 297 1173005 15072
<none> 1172708 15074
- region_Northeast 1 1627 1174335 15075
- BirthRate 1 2269 1174977 15076
- medIncome 1 3287 1175995 15078
- PercentMarried 1 15155 1187863 15103
- region_South 1 18757 1191465 15110
- region_West 1 25688 1198396 15124
- PctUnemployed16_Over 1 39994 1212702 15153
- PctHS25_Over 1 104760 1277467 15280
Step: AIC=15072.13
deathRate ~ medIncome + popEst2015 + PercentMarried + PctHS25_Over +
PctUnemployed16_Over + BirthRate + region_Northeast + region_South +
region_West
Df Sum of Sq RSS AIC
- popEst2015 1 251 1173189 15071
<none> 1172938 15072
- region_Northeast 1 1537 1174475 15073
+ PctPrivateCoverage 1 230 1172708 15074
- BirthRate 1 2125 1175064 15074
- medIncome 1 6658 1179596 15084
- PercentMarried 1 15328 1188266 15102
- region_South 1 23971 1196909 15119
- region_West 1 26930 1199868 15125
- PctUnemployed16_Over 1 48377 1221315 15169
- PctHS25_Over 1 104776 1277714 15279
Step: AIC=15070.65
deathRate ~ medIncome + PercentMarried + PctHS25_Over + PctUnemployed16_Over +
BirthRate + region_Northeast + region_South + region_West
Df Sum of Sq RSS AIC
<none> 1173189 15071
- region_Northeast 1 1623 1174812 15072
+ popEst2015 1 251 1172938 15072
+ PctPrivateCoverage 1 184 1173005 15072
- BirthRate 1 2108 1175297 15073
- medIncome 1 7272 1180460 15084
- PercentMarried 1 15080 1188269 15100
- region_South 1 24140 1197329 15118
- region_West 1 26985 1200174 15124
- PctUnemployed16_Over 1 48126 1221315 15167
- PctHS25_Over 1 109109 1282298 15285ajuste_AIC <-
lm(data = train,
formula =
deathRate ~ medIncome + PercentMarried + PctHS25_Over +
PctUnemployed16_Over + BirthRate + region_Northeast +
region_South + region_West)
ajuste_AIC |> summary()
Call:
lm(formula = deathRate ~ medIncome + PercentMarried + PctHS25_Over +
PctUnemployed16_Over + BirthRate + region_Northeast + region_South +
region_West, data = train)
Residuals:
Min 1Q Median 3Q Max
-103.450 -11.926 0.524 11.870 158.525
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.572e+02 6.065e+00 25.912 < 2e-16 ***
medIncome -2.012e-04 5.187e-05 -3.879 0.000108 ***
PercentMarried -4.713e-01 8.436e-02 -5.587 2.58e-08 ***
PctHS25_Over 1.238e+00 8.236e-02 15.027 < 2e-16 ***
PctUnemployed16_Over 1.705e+00 1.709e-01 9.980 < 2e-16 ***
BirthRate -4.851e-01 2.323e-01 -2.089 0.036829 *
region_Northeast -3.564e+00 1.945e+00 -1.833 0.066953 .
region_South 7.926e+00 1.121e+00 7.068 2.04e-12 ***
region_West -1.133e+01 1.517e+00 -7.473 1.09e-13 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 21.98 on 2428 degrees of freedom
Multiple R-squared: 0.3733, Adjusted R-squared: 0.3712
F-statistic: 180.8 on 8 and 2428 DF, p-value: < 2.2e-16ajuste_BIC <- MASS::stepAIC(ajuste_saturado, direction = "both",
k = log(nrow(train)))Start: AIC=15137.43
deathRate ~ medIncome + popEst2015 + PercentMarried + PctHS25_Over +
PctUnemployed16_Over + PctPrivateCoverage + BirthRate + region_Northeast +
region_South + region_West
Df Sum of Sq RSS AIC
- PctPrivateCoverage 1 230 1172938 15130
- popEst2015 1 297 1173005 15130
- region_Northeast 1 1627 1174335 15133
- BirthRate 1 2269 1174977 15134
- medIncome 1 3287 1175995 15136
<none> 1172708 15137
- PercentMarried 1 15155 1187863 15161
- region_South 1 18757 1191465 15168
- region_West 1 25688 1198396 15182
- PctUnemployed16_Over 1 39994 1212702 15211
- PctHS25_Over 1 104760 1277467 15338
Step: AIC=15130.11
deathRate ~ medIncome + popEst2015 + PercentMarried + PctHS25_Over +
PctUnemployed16_Over + BirthRate + region_Northeast + region_South +
region_West
Df Sum of Sq RSS AIC
- popEst2015 1 251 1173189 15123
- region_Northeast 1 1537 1174475 15126
- BirthRate 1 2125 1175064 15127
<none> 1172938 15130
- medIncome 1 6658 1179596 15136
+ PctPrivateCoverage 1 230 1172708 15137
- PercentMarried 1 15328 1188266 15154
- region_South 1 23971 1196909 15172
- region_West 1 26930 1199868 15178
- PctUnemployed16_Over 1 48377 1221315 15221
- PctHS25_Over 1 104776 1277714 15331
Step: AIC=15122.84
deathRate ~ medIncome + PercentMarried + PctHS25_Over + PctUnemployed16_Over +
BirthRate + region_Northeast + region_South + region_West
Df Sum of Sq RSS AIC
- region_Northeast 1 1623 1174812 15118
- BirthRate 1 2108 1175297 15119
<none> 1173189 15123
- medIncome 1 7272 1180460 15130
+ popEst2015 1 251 1172938 15130
+ PctPrivateCoverage 1 184 1173005 15130
- PercentMarried 1 15080 1188269 15146
- region_South 1 24140 1197329 15165
- region_West 1 26985 1200174 15170
- PctUnemployed16_Over 1 48126 1221315 15213
- PctHS25_Over 1 109109 1282298 15332
Step: AIC=15118.41
deathRate ~ medIncome + PercentMarried + PctHS25_Over + PctUnemployed16_Over +
BirthRate + region_South + region_West
Df Sum of Sq RSS AIC
- BirthRate 1 1742 1176554 15114
<none> 1174812 15118
+ region_Northeast 1 1623 1173189 15123
+ popEst2015 1 337 1174475 15126
+ PctPrivateCoverage 1 98 1174714 15126
- medIncome 1 9521 1184333 15130
- PercentMarried 1 13708 1188520 15139
- region_West 1 25366 1200178 15163
- region_South 1 31422 1206234 15175
- PctUnemployed16_Over 1 47178 1221990 15207
- PctHS25_Over 1 107638 1282450 15324
Step: AIC=15114.22
deathRate ~ medIncome + PercentMarried + PctHS25_Over + PctUnemployed16_Over +
region_South + region_West
Df Sum of Sq RSS AIC
<none> 1176554 15114
+ BirthRate 1 1742 1174812 15118
+ region_Northeast 1 1257 1175297 15119
+ popEst2015 1 308 1176246 15121
+ PctPrivateCoverage 1 26 1176528 15122
- medIncome 1 8811 1185365 15125
- PercentMarried 1 15444 1191998 15138
- region_West 1 25515 1202069 15159
- region_South 1 32512 1209067 15173
- PctUnemployed16_Over 1 47741 1224296 15203
- PctHS25_Over 1 109486 1286040 15323ajuste_BIC <-
lm(data = train,
formula =
deathRate ~ medIncome + PercentMarried + PctHS25_Over + PctUnemployed16_Over + region_South + region_West)
ajuste_BIC |> summary()
Call:
lm(formula = deathRate ~ medIncome + PercentMarried + PctHS25_Over +
PctUnemployed16_Over + region_South + region_West, data = train)
Residuals:
Min 1Q Median 3Q Max
-105.425 -12.091 0.501 11.766 159.447
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.541e+02 5.944e+00 25.933 < 2e-16 ***
medIncome -2.141e-04 5.020e-05 -4.266 2.07e-05 ***
PercentMarried -4.624e-01 8.187e-02 -5.648 1.82e-08 ***
PctHS25_Over 1.233e+00 8.202e-02 15.038 < 2e-16 ***
PctUnemployed16_Over 1.694e+00 1.706e-01 9.930 < 2e-16 ***
region_South 8.706e+00 1.062e+00 8.194 4.02e-16 ***
region_West -1.072e+01 1.476e+00 -7.259 5.21e-13 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 22 on 2430 degrees of freedom
Multiple R-squared: 0.3715, Adjusted R-squared: 0.37
F-statistic: 239.4 on 6 and 2430 DF, p-value: < 2.2e-16performance::compare_performance(ajuste_saturado, ajuste_AIC, ajuste_BIC)# Comparison of Model Performance Indices
Name | Model | AIC (weights) | AICc (weights) | BIC (weights) | R2 | R2 (adj.) | RMSE | Sigma
-------------------------------------------------------------------------------------------------------------------
ajuste_saturado | lm | 21991.6 (0.154) | 21991.7 (0.151) | 22061.1 (<.001) | 0.374 | 0.371 | 21.936 | 21.986
ajuste_AIC | lm | 21988.6 (0.690) | 21988.6 (0.691) | 22046.5 (0.013) | 0.373 | 0.371 | 21.941 | 21.982
ajuste_BIC | lm | 21991.5 (0.156) | 21991.6 (0.158) | 22037.9 (0.987) | 0.372 | 0.370 | 21.972 | 22.004Ajuste AIC: ha eliminado las variables
popEst2015yPctPrivateCoverage, empeorando ligeramente \(R^2\) pero incrementando ligeramente el \(R^2\) (menos variables pero más necesarias)Ajuste BIC: ha eliminado además las variables
BirthRateyregion_Northeast, empeorando ligeramente \(R^2\)
La diferencia de ambos es que el segundo es más “agresivo” penalizando. ¿La ventaja? Que asintóticamente (si \(n\) fuese incrementándose) tenemos garantizado que el BIC acaba eligiendo el modelo real que subyace en los datos (consistencia)
Al final tenemos 3 modelos con prácticamente la misma cantidad de información explicada, pero uno tiene 10 variables, otro 8 y otro 6 variables (con casi la misma calidad, pero con casi la mitad de predictoras) –> nos qedamos con el ajuste BIC para la diagnosis.
Ejercicio 6 (1.5 puntos)
Compara la calidad de los 3 modelos (saturado, BIC y AIC) en train. Quédate con uno (justifica) y comprueba si se cumplen las hipotesis necesarias (fase de diagnosis SOLO en ese elegido). ¿Se cumplen las hipótesis? Argumenta porqué, más allá de lo que valga luego la bondad de ajuste, es importante que se cumplan.
# Escribe el código que consideres
performance::check_model(ajuste_BIC)
En un primer vistazo visual vemos que es bastante probable que muchas hipótesis se cumplan o estén derca de complirse (teniendo en cuenta que no hemos tratado los outliers que pueden estar influenciando)
Linealidad
ggplot(tibble("fitted" = ajuste_BIC$fitted.values,
"residuals" = ajuste_BIC$residuals),
aes(x = fitted, y = residuals)) +
geom_point() +
geom_smooth(method = "lm") +
theme_minimal()`geom_smooth()` using formula = 'y ~ x'
A la izquierda y arriba se observan algunas observaciones que probablemente sean outliers peor el resto si parece presentar un patrón independiente entre residuos y valores ajustados, no viéndose ninguna tendencia o modelo posible a ajustar entre ellos, así que podríamos dar por válida la linealidad
Homocedasticidad
performance::check_heteroscedasticity(ajuste_BIC)Warning: Heteroscedasticity (non-constant error variance) detected (p < .001).ggplot(tibble("id" = 1:length(ajuste_BIC$residuals),
"residuals" = ajuste_BIC$residuals),
aes(x = id, y = residuals)) +
geom_point() +
geom_smooth(method = "lm") +
theme_minimal()`geom_smooth()` using formula = 'y ~ x'
A pesar de que el p-valor nos dice que se rechaz ala homocedasticidad, los errores no presentan ningún patrón ni aumentan su varianza, estando prácticamente todos dentro de una banda constante (salvo unos que probablemente sean outliers y deberíamos analizarlos de manera particular)
Normalidad
performance::check_normality(ajuste_BIC)Warning: Non-normality of residuals detected (p < .001).ggplot(tibble("residuals" = ajuste_BIC$residuals)) +
geom_qq(aes(sample = residuals)) +
geom_qq_line(aes(sample = residuals)) +
theme_minimal()
Normalidad no parece cumplirse. Recuerda que aquí influye:
- tamaño muestral
- outliers que hace que en los extremos de los quantiles la distribució se aleje
ggplot(tibble("residuals" = ajuste_BIC$residuals,
"normal" = rnorm(n = length(ajuste_BIC$residuals), mean = 0,
sd = sd(ajuste_BIC$residuals)))) +
geom_density(aes(x = residuals), color = "blue") +
geom_density(aes(x = normal), color = "red") +
theme_minimal()
Vemos como la distribución de los residuos es algo más apuntada que lo que correspondería a una distribución normal (probablemente por los outliers que contienen nuestros datos)
Correlación
performance::check_autocorrelation(ajuste_BIC)OK: Residuals appear to be independent and not autocorrelated (p = 0.612).acf(ajuste_BIC$residuals)
Vemos tanto numérica como visualmente uqe los residuos están incorrelados (con acf() debemos obtener lo contrario a lo que tendríamos al inicio de un análisis de series temporales: los residuos retardados sin correlacion entre ellos)
Predicciones simulación
performance::check_predictions(ajuste_BIC)Warning: Maximum value of original data is not included in the
replicated data.
Model may not capture the variation of the data.
Las simulaciones (bajo la hipotesis de que el modelo ajustado fuese el modelo real subyacente) están ligeramente por debajo del dato observado, probablemente producido por ese apuntamiento que se observaba en los residuales estimados.
Ejercicio 7 (0.75 puntos)
Asume que se cumplen las hipótesis (en caso de que no se cumpliesen) e interpreta todo lo que sepas sobre la segunda tabla de la salida de la regresión (la de inferencia de los parámetros). Elimina variables si su efecto no fuese significativo.
ajuste_BIC |> summary()
Call:
lm(formula = deathRate ~ medIncome + PercentMarried + PctHS25_Over +
PctUnemployed16_Over + region_South + region_West, data = train)
Residuals:
Min 1Q Median 3Q Max
-105.425 -12.091 0.501 11.766 159.447
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.541e+02 5.944e+00 25.933 < 2e-16 ***
medIncome -2.141e-04 5.020e-05 -4.266 2.07e-05 ***
PercentMarried -4.624e-01 8.187e-02 -5.648 1.82e-08 ***
PctHS25_Over 1.233e+00 8.202e-02 15.038 < 2e-16 ***
PctUnemployed16_Over 1.694e+00 1.706e-01 9.930 < 2e-16 ***
region_South 8.706e+00 1.062e+00 8.194 4.02e-16 ***
region_West -1.072e+01 1.476e+00 -7.259 5.21e-13 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 22 on 2430 degrees of freedom
Multiple R-squared: 0.3715, Adjusted R-squared: 0.37
F-statistic: 239.4 on 6 and 2430 DF, p-value: < 2.2e-16Asumiendo que las hipótesis son ciertas, interpretando los p-valores, vemos todas las covariables presentan un efecto lineal significativo así como a nivel global según el ANOVA realizado
ajuste_BIC |> anova()Analysis of Variance Table
Response: deathRate
Df Sum Sq Mean Sq F value Pr(>F)
medIncome 1 362519 362519 748.729 < 2.2e-16 ***
PercentMarried 1 26011 26011 53.722 3.128e-13 ***
PctHS25_Over 1 152438 152438 314.839 < 2.2e-16 ***
PctUnemployed16_Over 1 53142 53142 109.757 < 2.2e-16 ***
region_South 1 75831 75831 156.618 < 2.2e-16 ***
region_West 1 25515 25515 52.697 5.214e-13 ***
Residuals 2430 1176554 484
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Así lo confirmamos con el anova secuencial: en cada paso, se calcula la diferencia de SSR entre el modelo con las variables anteriroes y el modelo con todas más la variable correspondiente
Ejercicio 8 (1.5 puntos)
Por último, utiliza el dataset test para, aplicando los 3 modelos entrenados (saturado, AIC y BIC), predecir las observaciones de test (con cada uno de los modelos. Con las observaciones de test construye un nuevo dataset de 3 columnas: la variable objetivo, la predicción y el modelo usado (saturado, BIC, AIC). Calcula el \(R^2\) en el dataset de test para cada modelo. ¿Cuál funciona mejor en test? ¿Qué % de información explica cada modelo? ¿Cuánto vale su varianza residual?
# Escribe el código que consideres
pred <-
tibble("modelo" = "saturado",
"real" = test$deathRate,
"pred" = predict(ajuste_saturado, test)) |>
bind_rows(tibble("modelo" = "AIC",
"real" = test$deathRate,
"pred" = predict(ajuste_AIC, test))) |>
bind_rows(tibble("modelo" = "BIC",
"real" = test$deathRate,
"pred" = predict(ajuste_BIC, test)))
pred |>
summarise(SSR = sum((real - pred)^2),
SST = sum((real - mean(real))^2),
var_res = SSR / nrow(train),
R2 = 1 - SSR/SST, .by = modelo)# A tibble: 3 × 5
modelo SSR SST var_res R2
<chr> <dbl> <dbl> <dbl> <dbl>
1 saturado 323375. 473810. 133. 0.317
2 AIC 324694. 473810. 133. 0.315
3 BIC 328353. 473810. 135. 0.307Vemos como el modelo con un \(R^2\) más bajo es el ajuste BIC, pero con muy poca diferencia (algo que ya intuíamos ya que sucedía así en train). En los 3 casos los modelos explican aprox el 30% de información contenida en la variable objetivo